/-
Copyright (c) 2023 David Loeffler. All rights reserved.
Released under Apache 2.0 license as described in the file LICENSE.
Authors: David Loeffler
-/
module

public import Mathlib.Analysis.Fourier.AddCircle
public import Mathlib.MeasureTheory.Integral.Pi

/-!
# Multivariate Fourier series

In this file we define the Fourier series of an L² function on the `d`-dimensional unit circle, and
show that it converges to the function in the L² norm. We also prove uniform convergence of the
Fourier series if `f` is continuous and the sequence of its Fourier coefficients is summable.
-/

@[expose] public section

noncomputable section

open scoped BigOperators ComplexConjugate ENNReal

open Set Algebra Submodule MeasureTheory

-- some instances for unit circle

attribute [local instance] Real.fact_zero_lt_one

/-- In this file we normalise the measure on `ℝ / ℤ` to have total volume 1. -/
local instance : MeasureSpace UnitAddCircle := ⟨AddCircle.haarAddCircle⟩

/-- The measure on `ℝ / ℤ` is a Haar measure. -/
local instance : Measure.IsAddHaarMeasure (volume : Measure UnitAddCircle) :=
  inferInstanceAs (Measure.IsAddHaarMeasure AddCircle.haarAddCircle)

/-- The measure on `ℝ / ℤ` is a probability measure. -/
local instance : IsProbabilityMeasure (volume : Measure UnitAddCircle) :=
  inferInstanceAs (IsProbabilityMeasure AddCircle.haarAddCircle)

/-- The product of finitely many copies of the unit circle, indexed by `d`. -/
abbrev UnitAddTorus (d : Type*) := d → UnitAddCircle

namespace UnitAddTorus

variable {d : Type*} [Fintype d]

section Monomials

variable (n : d → ℤ)

/-- Exponential monomials in `d` variables. -/
def mFourier : C(UnitAddTorus d, ℂ) where
  toFun x := ∏ i : d, fourier (n i) (x i)
  continuous_toFun := by fun_prop

variable {n} {x : UnitAddTorus d}

lemma mFourier_neg : mFourier (-n) x = conj (mFourier n x) := by
  simp only [mFourier, Pi.neg_apply, fourier_neg, ContinuousMap.coe_mk, map_prod]

lemma mFourier_add {m : d → ℤ} : mFourier (m + n) x = mFourier m x * mFourier n x := by
  simp only [mFourier, Pi.add_apply, fourier_add, ContinuousMap.coe_mk, ← Finset.prod_mul_distrib]

lemma mFourier_zero : mFourier (0 : d → ℤ) = 1 := by
  ext x
  simp only [mFourier, Pi.zero_apply, fourier_zero, Finset.prod_const_one, ContinuousMap.coe_mk,
    ContinuousMap.one_apply]

lemma mFourier_norm : ‖mFourier n‖ = 1 := by
  apply le_antisymm
  · refine (ContinuousMap.norm_le _ zero_le_one).mpr fun i ↦ ?_
    simp only [mFourier, fourier_apply, ContinuousMap.coe_mk, norm_prod, Circle.norm_coe,
      Finset.prod_const_one, le_rfl]
  · refine (le_of_eq ?_).trans ((mFourier n).norm_coe_le_norm fun _ ↦ 0)
    simp only [mFourier, ContinuousMap.coe_mk, fourier_eval_zero, Finset.prod_const_one,
      CStarRing.norm_one]

lemma mFourier_single [DecidableEq d] (z : d → AddCircle (1 : ℝ)) (i : d) :
    mFourier (Pi.single i 1) z = fourier 1 (z i) := by
  simp_rw [mFourier, ContinuousMap.coe_mk]
  have := Finset.prod_mul_prod_compl {i} (fun j ↦ fourier ((Pi.single i (1 : ℤ) : d → ℤ) j) (z j))
  rw [Finset.prod_singleton, Finset.prod_congr rfl (fun j hj ↦ ?_)] at this
  · rw [← this, Finset.prod_const_one, mul_one, Pi.single_eq_same]
  · rw [Finset.mem_compl, Finset.mem_singleton] at hj
    simp only [Pi.single_eq_of_ne hj, fourier_zero]

end Monomials

section Algebra

/-- The star subalgebra of `C(UnitAddTorus d, ℂ)` generated by `mFourier n` for `n ∈ ℤᵈ`. -/
def mFourierSubalgebra (d : Type*) [Fintype d] : StarSubalgebra ℂ C(UnitAddTorus d, ℂ) where
  toSubalgebra := Algebra.adjoin ℂ (range mFourier)
  star_mem' := by
    change Algebra.adjoin ℂ (range mFourier) ≤ star (Algebra.adjoin ℂ (range mFourier))
    refine adjoin_le ?_
    rintro _ ⟨n, rfl⟩
    refine subset_adjoin ⟨-n, ?_⟩
    ext1 x
    simp only [mFourier_neg, starRingEnd_apply, ContinuousMap.star_apply]

/-- The star subalgebra of `C(UnitAddTorus d, ℂ)` generated by `mFourier n` for `n ∈ ℤᵈ` is in fact
the linear span of these functions. -/
theorem mFourierSubalgebra_coe :
    (mFourierSubalgebra d).toSubalgebra.toSubmodule = span ℂ (range mFourier) := by
  apply adjoin_eq_span_of_subset
  refine .trans (fun x ↦ Submonoid.closure_induction (fun _ ↦ id) ⟨0, ?_⟩ ?_) subset_span
  · ext z
    simp only [mFourier, Pi.zero_apply, fourier_zero, Finset.prod_const, one_pow,
      ContinuousMap.coe_mk, ContinuousMap.one_apply]
  · rintro _ _ _ _ ⟨m, rfl⟩ ⟨n, rfl⟩
    refine ⟨m + n, ?_⟩
    ext z
    simp only [mFourier, Pi.add_apply, fourier_apply, fourier_add', Finset.prod_mul_distrib,
      ContinuousMap.coe_mk, ContinuousMap.mul_apply]

/-- The subalgebra of `C(UnitAddTorus d, ℂ)` generated by `mFourier n` for `n ∈ ℤᵈ` separates
points. -/
theorem mFourierSubalgebra_separatesPoints : (mFourierSubalgebra d).SeparatesPoints := by
  classical
  intro x y hxy
  rw [Ne, funext_iff, not_forall] at hxy
  obtain ⟨i, hi⟩ := hxy
  refine ⟨_, ⟨mFourier (Pi.single i 1), subset_adjoin ⟨Pi.single i 1, rfl⟩, rfl⟩, ?_⟩
  dsimp only
  rw [mFourier_single, mFourier_single, fourier_one, fourier_one, Ne, Subtype.coe_inj]
  contrapose! hi
  exact AddCircle.injective_toCircle one_ne_zero hi

/-- The subalgebra of `C(UnitAddTorus d, ℂ)` generated by `mFourier n` for `n : d → ℤ` is dense. -/
theorem mFourierSubalgebra_closure_eq_top : (mFourierSubalgebra d).topologicalClosure = ⊤ :=
  ContinuousMap.starSubalgebra_topologicalClosure_eq_top_of_separatesPoints _
    mFourierSubalgebra_separatesPoints

/-- The linear span of the monomials `mFourier n` is dense in `C(UnitAddTorus d, ℂ)`. -/
theorem span_mFourier_closure_eq_top :
    (span ℂ (range <| mFourier (d := d))).topologicalClosure = ⊤ := by
  rw [← mFourierSubalgebra_coe]
  exact congr_arg (Subalgebra.toSubmodule <| StarSubalgebra.toSubalgebra ·)
    mFourierSubalgebra_closure_eq_top

end Algebra

section Lp

/-- The family of monomials `mFourier n`, parametrized by `n : ℤᵈ` and considered as
elements of the `Lp` space of functions `UnitAddTorus d → ℂ`. -/
abbrev mFourierLp (p : ℝ≥0∞) [Fact (1 ≤ p)] (n : d → ℤ) :
    Lp ℂ p (volume : Measure (UnitAddTorus d)) :=
  ContinuousMap.toLp (E := ℂ) p volume ℂ (mFourier n)

theorem coeFn_mFourierLp (p : ℝ≥0∞) [Fact (1 ≤ p)] (n : d → ℤ) :
    mFourierLp p n =ᵐ[volume] mFourier n :=
  ContinuousMap.coeFn_toLp volume (mFourier n)

-- shortcut instance: avoids `synthInstance.maxHeartbeats` issue in `span_mFourierLp_closure_eq_top`
instance {p : ℝ≥0∞} [Fact (1 ≤ p)] :
    ContinuousSMul ℂ (Lp ℂ p (volume : Measure (UnitAddTorus d))) :=
  inferInstance

/-- For each `1 ≤ p < ∞`, the linear span of the monomials `mFourier n` is dense in the `Lᵖ` space
of functions on `UnitAddTorus d`. -/
theorem span_mFourierLp_closure_eq_top {p : ℝ≥0∞} [Fact (1 ≤ p)] (hp : p ≠ ∞) :
    (span ℂ (range (@mFourierLp d _ p _))).topologicalClosure = ⊤ := by
  simpa only [map_span, ContinuousLinearMap.coe_coe, ← range_comp, Function.comp_def] using
    (ContinuousMap.toLp_denseRange ℂ volume ℂ hp).topologicalClosure_map_submodule
      span_mFourier_closure_eq_top

/-- The monomials `mFourierLp 2 n` are an orthonormal set in `L²`. -/
theorem orthonormal_mFourier : Orthonormal ℂ (mFourierLp (d := d) 2) := by
  rw [orthonormal_iff_ite]
  intro m n
  simp only [ContinuousMap.inner_toLp, ← mFourier_neg, ← mFourier_add]
  split_ifs with h
  · simpa only [h, add_neg_cancel, mFourier_zero, probReal_univ, one_smul] using
      integral_const (α := UnitAddTorus d) (μ := volume) (1 : ℂ)
  rw [mFourier, ContinuousMap.coe_mk, MeasureTheory.integral_fintype_prod_volume_eq_prod]
  obtain ⟨i, hi⟩ := Function.ne_iff.mp h
  apply Finset.prod_eq_zero (Finset.mem_univ i)
  simpa only [eq_false_intro hi, if_false, ContinuousMap.inner_toLp, ← fourier_neg,
    ← fourier_add] using (orthonormal_iff_ite.mp <| orthonormal_fourier) (m i) (n i)

end Lp

section fourierCoeff

variable {E : Type} [NormedAddCommGroup E] [NormedSpace ℂ E]

/-- The `n`-th Fourier coefficient of a function `UnitAddTorus d → E`, for `E` a complete normed
`ℂ`-vector space, defined as the integral over `UnitAddTorus d` of `mFourier (-n) t • f t`. -/
def mFourierCoeff (f : UnitAddTorus d → E) (n : d → ℤ) : E := ∫ t, mFourier (-n) t • f t

end fourierCoeff

section FourierL2

local notation "L²(" α ")" => Lp ℂ 2 (volume : Measure α)

/-- We define `mFourierBasis` to be a `ℤᵈ`-indexed Hilbert basis for the `L²` space of functions
on `UnitAddTorus d`, which by definition is an isometric isomorphism from `L²(UnitAddTorus d)`
to `ℓ²(ℤᵈ, ℂ)`. -/
def mFourierBasis : HilbertBasis (d → ℤ) ℂ L²(UnitAddTorus d) :=
  HilbertBasis.mk orthonormal_mFourier (span_mFourierLp_closure_eq_top (by simp)).ge

/-- The elements of the Hilbert basis `mFourierBasis` are the functions `mFourierLp 2`, i.e. the
monomials `mFourier n` on `UnitAddTorus d` considered as elements of `L²`. -/
@[simp]
theorem coe_mFourierBasis : ⇑(mFourierBasis (d := d)) = mFourierLp 2 := HilbertBasis.coe_mk _ _

/-- Under the isometric isomorphism `mFourierBasis` from `L²(UnitAddTorus d)` to `ℓ²(ℤᵈ, ℂ)`,
the `i`-th coefficient is `mFourierCoeff f i`. -/
theorem mFourierBasis_repr (f : L²(UnitAddTorus d)) (i : d → ℤ) :
    mFourierBasis.repr f i = mFourierCoeff f i := by
  trans ∫ t, conj (mFourierLp 2 i t) * f t
  · rw [mFourierBasis.repr_apply_apply f i, MeasureTheory.L2.inner_def, coe_mFourierBasis]
    simp only [RCLike.inner_apply, mul_comm]
  · apply integral_congr_ae
    filter_upwards [coeFn_mFourierLp 2 i] with _ ht
    rw [ht, ← mFourier_neg, smul_eq_mul]

/-- The Fourier series of an `L2` function `f` sums to `f` in the `L²` norm. -/
theorem hasSum_mFourier_series_L2 (f : L²(UnitAddTorus d)) :
    HasSum (fun i ↦ mFourierCoeff f i • mFourierLp 2 i) f := by
  simpa [← coe_mFourierBasis, mFourierBasis_repr] using mFourierBasis.hasSum_repr f

/-- **Parseval's identity** for inner products: for `L²` functions `f, g` on `UnitAddTorus d`, the
inner product of the Fourier coefficients of `f` and `g` is the inner product of `f` and `g`. -/
theorem hasSum_prod_mFourierCoeff (f g : L²(UnitAddTorus d)) :
    HasSum (fun i ↦ conj (mFourierCoeff f i) * (mFourierCoeff g i)) (∫ t, conj (f t) * g t) := by
  simp_rw [mul_comm (conj _)]
  refine HasSum.congr_fun (mFourierBasis.hasSum_inner_mul_inner f g) (fun n ↦ ?_)
  simp only [← mFourierBasis_repr, HilbertBasis.repr_apply_apply, inner_conj_symm,
    mul_comm (inner ℂ f _)]

/-- **Parseval's identity** for norms: for an `L²` function `f` on `UnitAddTorus d`, the sum of the
squared norms of the Fourier coefficients equals the `L²` norm of `f`. -/
theorem hasSum_sq_mFourierCoeff (f : L²(UnitAddTorus d)) :
    HasSum (fun i ↦ ‖mFourierCoeff f i‖ ^ 2) (∫ t, ‖f t‖ ^ 2) := by
  simpa only [← RCLike.inner_apply', inner_self_eq_norm_sq, ← integral_re
    (L2.integrable_inner f f)] using RCLike.hasSum_re ℂ (hasSum_prod_mFourierCoeff f f)

end FourierL2

section Convergence

variable (f : C(UnitAddTorus d, ℂ))

theorem mFourierCoeff_toLp (n : d → ℤ) :
    mFourierCoeff (f.toLp 2 volume ℂ) n = mFourierCoeff f n :=
  integral_congr_ae (ae_eq_rfl.mul <| f.coeFn_toAEEqFun _)

variable {f}

/-- If the sequence of Fourier coefficients of `f` is summable, then the Fourier series converges
uniformly to `f`. -/
theorem hasSum_mFourier_series_of_summable (h : Summable (mFourierCoeff f)) :
    HasSum (fun i ↦ mFourierCoeff f i • mFourier i) f := by
  have sum_L2 := hasSum_mFourier_series_L2 (ContinuousMap.toLp 2 volume ℂ f)
  simp only [mFourierCoeff_toLp] at sum_L2
  refine ContinuousMap.hasSum_of_hasSum_Lp (.of_norm ?_) sum_L2
  simpa only [norm_smul, mFourier_norm, mul_one] using h.norm

/-- If the sequence of Fourier coefficients of `f` is summable, then the Fourier series of `f`
converges everywhere pointwise to `f`. -/
theorem hasSum_mFourier_series_apply_of_summable (h : Summable (mFourierCoeff f))
    (x : UnitAddTorus d) : HasSum (fun i ↦ mFourierCoeff f i • mFourier i x) (f x) := by
  simpa only [map_smul] using (ContinuousMap.evalCLM ℂ x).hasSum
    (hasSum_mFourier_series_of_summable h)

end Convergence

end UnitAddTorus
